/*
 * fft.cc
 *
 *  Created on: May 2, 2011
 *      Author: ilewis
 */
#include <complex>
#include <math.h>
#include <android/log.h>
#include <arm_neon.h>
using namespace std;

// A completely naive slow DFT. The only reason this exists is to check results from
// our more interesting fast implementations.
void NaiveFT(float* inR, float* inI, float* outR, float* outI, int log2_length) {
	int length = 1 << log2_length;
	complex<float>* out = new complex<float>[length];

	for (int n = 0; n < length; ++n) {
		out[n] = complex<float>(0,0);
		for (int k = 0; k < length; ++k) {
			out[n] += complex<float>(inR[k], inI[k])
			             * complex<float>(
			            		 cos((float)n*(float)k*2.0f*M_PI/length),
			            		 sin((float)n*(float)k*2.0f*M_PI/length));
		}
	}

}

int bitreverse(int in, int len) {
	// naive bitreverse
	int out = 0;
	for (int i = 0; i < len; ++i) {
		out <<= 1;
		out |= in & 1;
		in >>= 1;
	}
	return out;
}

void ComplexMult(float aR, float aI, float bR, float bI, float* oR, float* oI) {
	*oR = aR * bR - aI * bI;
	*oI = aR * bI + bR * aI;
}

// C++ model of FFT algo
void FFT_cpp(std::complex<float>* in, std::complex<float>* out, int log2_length) {
	int len = 1 << log2_length;

	for (int i = 0; i < len; ++i) {
		out[i] = in[i];
	}

	for (int n = 1; n <= log2_length; ++n) { 	// n = current stage
		int partitionSize = len / (1 << n); // size of each contiguous partition
		int nPartitions = len / partitionSize;

		__android_log_print(ANDROID_LOG_INFO, "DSP", "STAGE %d: partition size is %d", n, partitionSize);

		for (int j = 0; j < nPartitions / 2; ++j) {			// j = current partition
			__android_log_print(ANDROID_LOG_INFO, "DSP", " PART %d", j);
			for (int k = 0; k < partitionSize; ++k) {
				int indexA = j * 2 * partitionSize + k;
				int indexB = (j * 2 + 1) * partitionSize +k;
				complex<float> a = out[indexA];
				complex<float> b = out[indexB];


				float step = (float)k / ((float)partitionSize * 2.0f);
				float angle = 2.0f * M_PI * step;
				out[indexA] = a + b;
				out[indexB] = (a - b) * complex<float>(cos(angle), sin(angle));
				__android_log_print(ANDROID_LOG_INFO, "DSP", "   %d, %d: %f\n", indexA, indexB, step);
			}
		}
	}
	for (int i = 0; i < len / 2; ++i) {
		__android_log_print(ANDROID_LOG_INFO, "DSP", "%d => %d", i, bitreverse(i,log2_length));
		std::swap(out[i], out[bitreverse(i,log2_length)]);
	}
}

// BEGIN NEON implementation

